/**
 * This class describes lattice that we use for our Ising model
 */
package simulation;

import java.util.Random;

import org.apache.commons.math.complex.Complex;

import common.Neighbors;
import common.TupleXY;
import common.TupleXYZ;

/**
 * @author vladimir
 * 
 */
public class Lattice {
	int nSites = -1;

	public Interval[][] lattice2 = null;
	public Interval[] lattice1 = null;

	int[][][][] neighbors2 = null;
	int[][] neighbors1 = null;

	public int coordNumber = 0;
	int dimension = -1;
	double J = 0;
	double B = -1;
	int L = -1;
	double T = -1;
	double beta = -1;
	Random rand = new Random();

	TupleXY[] sublattice1 = null;
	TupleXY[] sublattice2 = null;
	TupleXY[] sublattice3 = null;

	public Lattice(int[][] tneighbors, double J, int L, double B, double T) {
		this.coordNumber = tneighbors.length;// number of the neighbors
		this.dimension = tneighbors[0].length;// dimensionality of the lattice
		this.nSites = (int) Math.pow(L, this.dimension);// total number of sites
														// in spatial direction
		this.L = L;// size of the lattice
		this.J = J;// colupling constant in spatial direction
		this.T = T;// temperature
		this.beta = 1 / T;
		this.B = B;

		if (this.dimension == 2) {
			Neighbors NB = new Neighbors(tneighbors, this.L, this.L);
			this.neighbors2 = NB.neighbors2;
			SubLattice Sub = new SubLattice(this.J, tneighbors, this.L, this.L);
			this.sublattice1 = Sub.sub1;
			this.sublattice2 = Sub.sub2;
			this.sublattice3 = Sub.sub3;
		} else if (this.dimension == 1) {
			Neighbors NB = new Neighbors(tneighbors, this.L);
			this.neighbors1 = NB.neighbors1;
		}

		if (dimension == 2) {
			this.lattice2 = new Interval[L][L];
			// create lattice array
			for (int x = 0; x < L; x++) {
				for (int y = 0; y < L; y++) {
					// fill lattice
					lattice2[x][y] = new Interval(B, T);
				}
			}
		} else if (dimension == 1) {
			this.lattice1 = new Interval[L];
			for (int x = 0; x < this.L; x++) {
				// fill lattice
				lattice1[x] = new Interval(B, T);
			}
		}
	}

	/*
	 * Calculates Energy of the Lattice
	 */
	public double Energy() {
		double result = 0;

		if (this.dimension == 2) {
			for (int x = 0; x < this.L; x++) {
				for (int y = 0; y < this.L; y++) {
					for (TupleXYZ cluster : lattice2[x][y].clusters) {
						result -= cluster.z
								* sumNeighbors(x, y, cluster.x, cluster.y);
					}
				}
			}
		} else if (this.dimension == 1) {
			for (int x = 0; x < this.L; x++) {
				for (TupleXYZ cluster : lattice1[x].clusters) {
					result -= cluster.z * sumNeighbors(x, cluster.x, cluster.y);
				}
			}
		}
		if (this.B != 0) {
			return result * this.J / 2 / this.nSites / this.beta - nCuts()
					/ this.beta;
		} else {
			return result * this.J / 2 / this.nSites / this.beta;
		}
	}

	/*
	 * Calculates Magnetisation of the lattice
	 */
	public double M(int power) {
		double result = 0;
		double[] percent_case = { 0, 0.25, 0.5, 0.75 };
		for (double percent : percent_case) {
			result += Math.pow(sum_over_slice(percent), power);
		}
		return result / percent_case.length / Math.pow(this.nSites, power);
	}

	public int sum_over_slice(double percent) {
		int result = 0;
		if (this.dimension == 2) {
			for (int x = 0; x < this.L; x++) {
				for (int y = 0; y < this.L; y++) {
					int index = this.lattice2[x][y].index(percent * this.beta);
					result += this.lattice2[x][y].clusters.get(index).z;
				}
			}
		} else if (this.dimension == 1) {
			for (int x = 0; x < this.L; x++) {
				int index = this.lattice1[x].index(percent * this.beta);
				result += this.lattice1[x].clusters.get(index).z;
			}
		}
		return result;
	}

	public double M_total() {
		double result = 0;
		if (this.dimension == 2) {
			for (int x = 0; x < this.L; x++) {
				for (int y = 0; y < this.L; y++) {
					result += this.lattice2[x][y].value();
				}
			}
		} else if (this.dimension == 1) {
			for (int x = 0; x < this.L; x++) {
				result += this.lattice1[x].value();
			}
		}
		return result / this.nSites / this.beta;
	}

	// calculates sum of neighbors within range not normalized. 2d case
	public double sumNeighbors(int siteX, int siteY, double lower, double upper) {
		double result = 0;
		for (int tn = 0; tn < this.coordNumber; tn++) {
			int neighborX = this.neighbors2[siteX][siteY][tn][0];
			int neighborY = this.neighbors2[siteX][siteY][tn][1];
			result += this.lattice2[neighborX][neighborY].value(lower, upper);
		}
		return result;
	}

	// calculates sum of neighbors within range not normalized. 1d case
	public double sumNeighbors(int siteX, double lower, double upper) {
		double result = 0;
		for (int tn = 0; tn < this.coordNumber; tn++) {
			int neighborX = this.neighbors1[siteX][tn];
			result += this.lattice1[neighborX].value(lower, upper);
		}
		return result;
	}

	// updates cluster
	public void update() {
		if (this.dimension == 2) {
			for (int x = 0; x < this.L; x++) {
				for (int y = 0; y < this.L; y++) {
					TupleXYZ newCluster = this.lattice2[x][y].newCluster();
					double dE = DeltaE(x, y, newCluster);
					if (B != 0) {
						if (rand.nextFloat() < Math.exp(-dE)) {
							this.lattice2[x][y].set(newCluster);
						}
					} else {
						if (rand.nextFloat() < 1 / (1 + Math.exp(dE))) {
							this.lattice2[x][y].set(newCluster);
						}
					}
				}
			}
		} else if (this.dimension == 1) {
			for (int x = 0; x < this.L; x++) {
				int tx = rand.nextInt(this.L);
				if (this.B != 0) {
					TupleXYZ newCluster = this.lattice1[tx].newCluster();
					double dE = DeltaE(tx, newCluster);
					if (rand.nextFloat() < Math.exp(-dE)) {
						this.lattice1[tx].set(newCluster);
					}
				} else {
					TupleXYZ newCluster = new TupleXYZ(0, 0, 0);
					double dE = DeltaE(tx, newCluster);
					if (rand.nextFloat() < 1 / (1 + Math.exp(dE))) {
						this.lattice1[tx].set(newCluster);
					}
				}
			}
		}

	}

	// calculates change in Energy between old and new configuration
	public double DeltaE(int x, int y, TupleXYZ newCluster) {
		return 2 * J * this.lattice2[x][y].clusters.get(newCluster.z).z
				* sumNeighbors(x, y, newCluster.x, newCluster.y);
	}

	// calculates change in Energy between old and new configuration
	public double DeltaE(int x, TupleXYZ newCluster) {
		int clustValue = this.lattice1[x].clusters.get(newCluster.z).z;
		if (this.B != 0) {
			return 2 * J * clustValue
					* sumNeighbors(x, newCluster.x, newCluster.y);
		} else {
			return 2 * J * clustValue * sumNeighbors(x, 0, 0);
		}
	}

	// returns average number of cuts
	public double nCuts() {
		double result = 0;
		if (this.dimension == 2) {
			for (int x = 0; x < this.L; x++) {
				for (int y = 0; y < this.L; y++) {
					result += this.lattice2[x][y].nCuts();
				}
			}
		} else if (this.dimension == 1) {
			for (int x = 0; x < this.L; x++) {
				result += this.lattice1[x].nCuts();
			}
		}
		return result / this.nSites;
	}

	// magnetisation on sublattice (index - index of the sublattice)
	public double M_sublattice(int sub_index) {
		if (this.J > 0) {
			return 0;
		}
		double result = 0;
		TupleXY[] temp = this.sublattice1;
		if (sub_index == 0) {
			temp = this.sublattice1;
		} else if (sub_index == 1) {
			temp = this.sublattice2;
		} else if (sub_index == 2) {
			temp = this.sublattice3;
		} else {
			System.err.println("wrong sub index : " + sub_index);
			System.exit(0);
		}
		for (TupleXY sub : temp) {
			result += this.lattice2[sub.x][sub.y].clusters.get(0).z;
		}

		return result / temp.length;
	}

	// total magnetisation on sublattice (index - index of the sublattice)
	public double M_sublattice_total(int sub_index) {
		if (J > 0) {
			return 0;
		}
		double result = 0;
		TupleXY[] temp = this.sublattice1;
		if (sub_index == 0) {
			temp = this.sublattice1;
		} else if (sub_index == 1) {
			temp = this.sublattice2;
		} else if (sub_index == 2) {
			temp = this.sublattice3;
		} else {
			System.err.println("wrong sub index : " + sub_index);
			System.exit(0);
		}

		for (TupleXY sub : temp) {
			result += this.lattice2[sub.x][sub.y].value();
		}

		return result / this.beta / temp.length;

	}

	// order parameter
	public double orderParameter() {
		double M0 = M_sublattice(0);
		double M1 = M_sublattice(1);
		double M2 = M_sublattice(2);

		return Math.sqrt(Math.abs(M0 * M0 + M1 * M1 + M2 * M2 - M0 * M1 - M0
				* M2 - M1 * M2) / 3);

	}

	// order parameter total
	public double orderParameter_total() {
		double M0 = M_sublattice_total(0);
		double M1 = M_sublattice_total(1);
		double M2 = M_sublattice_total(2);

		return Math.sqrt(Math.abs(M0 * M0 + M1 * M1 + M2 * M2 - M0 * M1 - M0
				* M2 - M1 * M2) / 3);
	}

	// Cos(6 * teta)
	public double cos6teta() {
		double M0 = M_sublattice(0);
		double M1 = M_sublattice(1);
		double M2 = M_sublattice(2);

		Complex phase1 = new Complex(-0.5, Math.sqrt(3) / 2);
		Complex phase2 = new Complex(-0.5, -Math.sqrt(3) / 2);

		Complex term0 = new Complex(M0, 0);
		Complex term1 = phase1.multiply(M1);
		Complex term2 = phase2.multiply(M2);

		Complex complexOrder = (term0.add(term1).add(term2))
				.divide(new Complex(Math.sqrt(3), 0));
		double argument = complexOrder.getArgument();
		if (Double.isNaN(argument)) {
			argument = rand.nextFloat() * 2 * Math.PI;
		}
		return Math.cos(6 * argument);
	}

	// Cos(6 * teta) total
	public double cos6teta_total() {
		double M0 = M_sublattice_total(0);
		double M1 = M_sublattice_total(1);
		double M2 = M_sublattice_total(2);
		Complex phase1 = new Complex(-0.5, Math.sqrt(3) / 2);
		Complex phase2 = new Complex(-0.5, -Math.sqrt(3) / 2);

		Complex term0 = new Complex(M0, 0);
		Complex term1 = phase1.multiply(M1);
		Complex term2 = phase2.multiply(M2);

		Complex complexOrder = (term0.add(term1).add(term2))
				.divide(new Complex(Math.sqrt(3), 0));
		double argument = complexOrder.getArgument();
		if (Double.isNaN(argument)) {
			argument = rand.nextFloat() * 2 * Math.PI;
		}
		return Math.cos(6 * argument);
	}

}
